function [MTF2, acut] =texture_mtf(varargin)
%[mtf, acutance] = texture_MTF  
% Image texture evaluation:  Texture Modulation Transfer Function (MTF)
% and acutance computed from deadleaves test target images. The method is
% based on computing a noise-power spectrum for the 'dead leaves' area of
% test images. This is done for an input (e.g. test target) image and an
% output image, e.g. after image processing. The texture MTF is computed 
% as the ratio square-root of the input and output spectra. The acutance is
% a summary measure, which provides a visually weighted summary measure of 
% the image texture capture(or retained) by the system under test.
%
%This implementation includs the two-dimensional trend removal described
%in the reference below. It also performs the image noise subtraction
%described by McElvain, et al.(2010).
%
%For background and details on this method, please see,
% P.D.Burns, Refined Measurement of Digital Image Texture Loss, Proc.
% SPIE vol. 8653, Image Quality and System Performance X, 86530H (2013) 
%
%Needed:
%Two corresponding input and output images files for the system being
%tested. Also,
% deadleavesNPS, radialNPS, texture_spec, rgb2lum, detrend2
%
% Peter D. Burns, pdburns@ieee.org 5 Nov. 2015

[f1,p1] = uigetimage('Select first (input) image file of texture target');
[Respec_sig1, freq1] = texture_spec([p1,f1], 1); 

[f2,p2] = uigetimage('Select second (output) image file of texture target');
[Respec_sig2, freq2] = texture_spec([p2,f2], 1);

[nn1]=length(Respec_sig1);
[nn2]=length(Respec_sig2);
if nn2~=nn1
    Rspec_sig2 = interp1(freq2, Respec1_sig2, freq1, 'spline');
end

% Compute the texture MTF as ssquare root of ratio of the spectra
MTF = (Respec_sig2./Respec_sig1);
MTF(1)=1;
MTF = sqrt(abs(MTF));

% Low-frequency normalization (scaling) at a frequency = 0.02 cy/pixel
ii = find(abs(freq1-0.02) == min(abs(freq1-0.02)));
MTF2 = MTF/MTF(ii);
MTF2(1:ii) = 1;

% Acutance Measurement
dx = 2.54/100;  % 100 ppi display
vd = 60; %cm     % viewing distance
csf = csf1(freq1, vd, dx);
acut = sum(MTF2.*csf)/sum(csf); 

set(0, 'DefaultTextInterpreter', 'tex')
pos = centerfig(10,4,0.8);
h = figure('Position',pos);
subplot(1,2,1)
semilogy(freq1, Respec_sig1,'LineWidth',1.5), hold on
semilogy(freq1, Respec_sig2,'r--','LineWidth',1.5)
hold off
mmax = 1.05*max(Respec_sig1(:));
mmin = mmax*1e-6;
axis([0 .69 mmin mmax])
xlabel('Frquency, cy/pixel')
ylabel('Power Spectrum')
legend('Input','Output')

subplot(1,2,2)
plot(freq1, MTF2,'LineWidth',1.5);
xlabel('Frquency, cy/pixel')
ylabel('MTF_{txt}')
title(['Acutance: ',num2str(acut,3)])
axis([0 .69 0 1.5])

